home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / AA_51.ZIP / SIDRLT.C < prev    next >
C/C++ Source or Header  |  1993-02-09  |  2KB  |  74 lines
  1.  
  2. /* Local Apparent Sidereal Time with equation of the equinoxes
  3.  * AA page B6
  4.  *
  5.  * Caution. At epoch J2000.0, the 16 decimal precision
  6.  * of IEEE double precision numbers
  7.  * limits time resolution measured by Julian date
  8.  * to approximately 24 microseconds.
  9.  */
  10.  
  11.  
  12. extern double J2000, TDT, RTD, nutl, coseps;
  13. double floor();
  14. int nutlo(), epsiln();
  15.  
  16. /* program returns sidereal seconds since sidereal midnight */
  17. double sidrlt( j, tlong )
  18. double j;    /* Julian date and fraction */
  19. double tlong;    /* East longitude of observer, degrees */
  20. {
  21. double jd0;     /* Julian day at midnight Universal Time */
  22. double secs;   /* Time of day, UT seconds since UT midnight */
  23. double eqeq, gmst, jd, T0, T, msday;
  24. /*long il;*/
  25.  
  26.   /* Julian day at given UT */
  27. jd = j;
  28. jd0 = floor(jd);
  29. secs = j - jd0;
  30. if( secs < 0.5 )
  31.     {
  32.     jd0 -= 0.5;
  33.     secs += 0.5;
  34.     }
  35. else
  36.     {
  37.     jd0 += 0.5;
  38.     secs -= 0.5;
  39.     }
  40. secs *= 86400.0;
  41.  
  42.   /* Julian centuries from standard epoch J2000.0 */
  43. T = (jd - J2000)/36525.0;
  44.   /* Same but at 0h Universal Time of date */
  45. T0 = (jd0 - J2000)/36525.0;
  46.  
  47. /* The equation of the equinoxes is the nutation in longitude
  48.  * times the cosine of the obliquity of the ecliptic.
  49.  * We already have routines for these.  Using the nutation formula
  50.  * given by Meeus, the peak error in 1986 is 5 milliseconds.
  51.  */
  52. nutlo(TDT);
  53. epsiln(TDT);
  54. /* nutl is in radians; convert to seconds of time
  55.  * at 240 seconds per degree
  56.  */
  57. eqeq = 240.0 * RTD * nutl * coseps;
  58.   /* Greenwich Mean Sidereal Time at 0h UT of date */
  59. gmst = (( -6.2e-6*T0 + 9.3104e-2)*T0 + 8640184.812866)*T0 + 24110.54841;
  60. /* mean solar days per sidereal day at date T0, = 1.00273790934 in 1986 */
  61. msday = 1.0 + ((-1.86e-5*T0 + 0.186208)*T0 + 8640184.812866)/(86400.*36525.);
  62.   /* Local apparent sidereal time at given UT */
  63. gmst = gmst + msday*secs + eqeq + 240.0*tlong;
  64.   /* Sidereal seconds modulo 1 sidereal day */
  65. gmst = gmst - 86400.0 * floor( gmst/86400.0 );
  66. /*
  67.  * il = gmst/86400.0;
  68.  * gmst = gmst - 86400.0 * il;
  69.  * if( gmst < 0.0 )
  70.  *    gmst += 86400.0;
  71.  */
  72. return( gmst );
  73. }
  74.